############################################################################################################################################
Downstream data analysis using MicrobiomeAnalyst web server or MicrobiomeAnalyst R package
RESISPART: Bioinformatics Workshop Tutorial, 2019, Brazil
Achal Dhariwal, Tsute Chen (George)
############################################################################################################################################

#######################################################################################################
#1. MicrobiomeAnalyst web server (Dhariwal et al., 2017)
MicrobiomeAnalyst is freely available at http://www.microbiomeanalyst.ca
Article can be accessed from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5570177/
                   OR
#2. MicrobiomeAnalystR: companion R package (under construction)
Link: https://github.com/xia-lab/MicrobiomeAnalystR
###############################################
#Step 1. Install R packages
#CRAN packages
install.packages("randomForest")
install.packages("viridis")
install.packages("RColorBrewer")
install.packages("pheatmap")
install.packages("Cairo")
install.packages("data.table")
install.packages("pheatmap")
install.packages("reshape")
install.packages("splitstackshape")
#Bioconductor packages
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("metagenomeSeq")
BiocManager::install("phyloseq") #if not already installed previously

#Step2: Open the R script name (MA_backendR.R){Open Script in R GUI}
#Step3: Run all the code blocks {Ctrl+A (Select All) then Ctrl+R (Run)}
###############################################
#######################################################################################################

############Downstream analysis########################################################################

###############################################
#initaite object
mbSet<-Init.mbSetObj()
#set the analysis module
mbSet<-SetModuleType(mbSet, "mdp")

####Data Input and Processing ################
#Step 1: reading files (output from dada2)
mbSet<-ReadSampleTable(mbSet, "meta.txt") #metadata
mbSet<-Read16STaxaTable(mbSet, "tax_table.txt") #taxonomy table
mbSet<-Read16SAbundData(mbSet, "otu_table.txt","text","Others/Not_specific","T") #ASV/OTU table
#Step 2: Data Sanity checking
mbSet<-SanityCheckData(mbSet, "text")
mbSet<-PlotLibSizeView(mbSet, "norm_libsizes_0","png")
#intermediate step (Create phyloseq object)
mbSet<-CreatePhyloseqObj(mbSet, "text","Others/Not_specific","F")
#Step 3: Data filteration
mbSet<-ApplyAbundanceFilter(mbSet, "prevalence", 4, 0.2)
mbSet<-ApplyVarianceFilter(mbSet, "iqr", 0.1)
#Step 4: Data Normalization
mbSet<-PerformNormalization(mbSet, "none", "colsum", "none")
###############################################

#1. Abundance bar plot (at Phylum, Genus, Species)
mbSet<-PlotTaxaAundanceBar(mbSet, "taxa_alpha_1","Phylum","Concentration", "Concentration", "barnorm",10, "grad","sum","png")
mbSet<-PlotTaxaAundanceBar(mbSet, "taxa_alpha_2","Genus","Concentration", "Concentration", "barnorm",10, "grad","sum","png")
mbSet<-PlotTaxaAundanceBar(mbSet, "taxa_alpha_3","Species","Concentration", "Concentration", "barnorm",10, "grad","sum","png")
###############################################

#2. Alpha diversity
mbSet<-PlotAlphaData(mbSet, "orig","alpha_diver_3","Chao1","Concentration","OTU", "default", "png")
mbSet<-PlotAlphaBoxData(mbSet, "alpha_diverbox_3","Chao1","Concentration","default", "png")
###############################################

#3. Beta Diversity (2D PCoA)
mbSet<-PlotBetaDiversity(mbSet, "beta_diver_0","PCoA","bray","expfac","Concentration","none","OTU","","Chao1", "yes", "png", 72, "default")
###############################################

#4. Rarefaction curves
mbSet<-PlotRarefactionCurve(mbSet, "filt","Concentration","Concentration","Concentration","5","rarefaction_curve_0","png")
###############################################

#5. Heatmap (Clustering analysis)
mbSet<-PlotHeatmap(mbSet, "heatmap_1","euclidean","ward.D","bwm","Concentration","Species","overview","F", "png","T","F")
###############################################

#6. Differential abundance testing (metagenomeSeq: R package)
mbSet<-PerformMetagenomeSeqAnal(mbSet, "Concentration",0.05,"NA","Species","zigfit")
###############################################

7. Random forest classification
mbSet<-RF.Anal(mbSet, 500,7,1,"Concentration","Genus")
mbSet<-PlotRF.Classify(mbSet, 15, "rf_cls_2","png", width=NA)
mbSet<-PlotRF.VIP(mbSet, 15, "rf_imp_2","png", width=NA)
###############################################
#Note: A lot of more functions are available for users to try and analyze thier dataset
#######################################################################################################
